#include "AnaDef.h"
#include "style.h"
#include "XGLUtils.h"

const Double_t EPSILON = 1e-10;

void downscaleHist(TH1D *hh, const Double_t fac)
{
  //take into account stat error of fac

  for(Int_t ii=1; ii<=hh->GetNbinsX(); ii++){
    const Double_t yy = hh->GetBinContent(ii);
    if(yy<EPSILON)
      continue;

    const Double_t newy = yy/fac;
    
    const Double_t newE = newy* TMath::Sqrt(1/yy+1/fac);

    hh->SetBinContent(ii, newy); hh->SetBinError(ii, newE);
  }
}

void normHist(TH1D* hh)
{
  for(Int_t ii=1; ii<=hh->GetNbinsX(); ii++){
    const Double_t cc = hh->GetBinContent(ii);
    const Double_t err = hh->GetBinError(ii);
    const Double_t bw = hh->GetBinWidth(ii);

    hh->SetBinContent(ii, cc/bw);
    hh->SetBinError(ii, err/bw);
  }
}

void getRawYield(TList *ll,const Int_t ktpc)
{
  TH1D * hjetpt= (TH1D*) ll->FindObject("fh1FFJetPtRecCutsInc");
  if(!hjetpt){
    printf("fh1FFJetPtRecCutsInc null!!\n");
    exit(1);
  }

  TH2D * htrkpt =(TH2D*) ll->FindObject("fh2FFTrackPtRecCutsInc");
  if(!htrkpt){
    printf("fh2FFTrackPtRecCutsInc null!!\n");
    exit(1);
  }

  TH2D * hz =(TH2D*) ll->FindObject("fh2FFZRecCutsInc");
  if(!hz){
    printf("fh2FFZRecCutsInc null!!\n");
    exit(1);
  }

  TH2D * hxi =(TH2D*) ll->FindObject("fh2FFXiRecCutsInc");
  if(!hxi){
    printf("fh2FFXiRecCutsInc null!!\n");
    exit(1);
  }

  hjetpt->RebinX(5);
  htrkpt->RebinX(5);
  hz->RebinX(5);
  hxi->RebinX(5);

  const Double_t jpt0[]={5,10,15};
  const Double_t jpt1[]={10,15,20};
  const Int_t jptb[]={1,2,3};
  Double_t njet[]={-999,-999,-999};

  for(Int_t ii = 0; ii < (Int_t)(sizeof(jptb)/sizeof(Int_t)); ii++){
    const Int_t ibin = jptb[ii];

    TList * lout=new TList;
    TH1D * htmp = 0x0;
    TGraphAsymmErrors * grtmp = 0x0;

    if( TMath::Abs(hjetpt->GetXaxis()->GetBinLowEdge(ibin)-jpt0[ii]) >EPSILON || 
        TMath::Abs(hjetpt->GetXaxis()->GetBinUpEdge(ibin) -jpt1[ii]) >EPSILON ){
      printf("bad jet bin in hjetpt !! ii %d ibin %d jpt0 %e jpt1 %e low %e up %e\n", ii, ibin, jpt0[ii], jpt1[ii], hjetpt->GetXaxis()->GetBinLowEdge(ibin), hjetpt->GetXaxis()->GetBinUpEdge(ibin)); exit(1);
    }
    njet[ii] = hjetpt->GetBinContent(ibin);

    //==================================

    if( TMath::Abs(htrkpt->GetXaxis()->GetBinLowEdge(ibin)-jpt0[ii]) >EPSILON || 
        TMath::Abs(htrkpt->GetXaxis()->GetBinUpEdge(ibin) -jpt1[ii]) >EPSILON ){
      printf("bad jet bin in htrkpt !! ii %d ibin %d jpt0 %e jpt1 %e low %e up %e\n", ii, ibin, jpt0[ii], jpt1[ii], htrkpt->GetXaxis()->GetBinLowEdge(ibin), htrkpt->GetXaxis()->GetBinUpEdge(ibin)); exit(1);
    }
    htmp = htrkpt->ProjectionY(Form("htrkpt%d",ii), ibin, ibin); lout->Add(htmp); 
    //if(htmp->GetSumw2()==0) {htmp->Sumw2();} htmp->Scale(1./njet[ii]); 
    downscaleHist(htmp, njet[ii]);
    style::ToNaturalScale(htmp); 
    //normHist after toNatural
    normHist(htmp);
    htmp->GetYaxis()->SetTitle("dN/dpt/Njet");

    //TGraphAsymmErrors* HistToGraph(const TString grname, const TH1* hh, const Double_t thres = 0, const TH1* herr = 0x0, const Double_t xmin = -1e10, const Double_t xmax = 1e10, const Bool_t klogx = kFALSE)
    grtmp = XGLUtils::HistToGraph("grid1", htmp, 0, 0x0, -1e10, 1e10, kTRUE);
    lout->Add(grtmp);


    if( TMath::Abs(hz->GetXaxis()->GetBinLowEdge(ibin)-jpt0[ii]) >EPSILON || 
        TMath::Abs(hz->GetXaxis()->GetBinUpEdge(ibin) -jpt1[ii]) >EPSILON ){
      printf("bad jet bin in hz !! ii %d ibin %d jpt0 %e jpt1 %e low %e up %e\n", ii, ibin, jpt0[ii], jpt1[ii], hz->GetXaxis()->GetBinLowEdge(ibin), hz->GetXaxis()->GetBinUpEdge(ibin)); exit(1);
    }
    htmp = hz->ProjectionY(Form("hz%d",ii), ibin, ibin); lout->Add(htmp); 
    //if(htmp->GetSumw2()==0) {htmp->Sumw2();} htmp->Scale(1./njet[ii]); 
    downscaleHist(htmp, njet[ii]);
    normHist(htmp);
    htmp->GetYaxis()->SetTitle("dN/dz/Njet");
    grtmp = XGLUtils::HistToGraph("grid2", htmp); lout->Add(grtmp);

    if( TMath::Abs(hxi->GetXaxis()->GetBinLowEdge(ibin)-jpt0[ii]) >EPSILON || 
        TMath::Abs(hxi->GetXaxis()->GetBinUpEdge(ibin) -jpt1[ii]) >EPSILON ){
      printf("bad jet bin in hxi !! ii %d ibin %d jpt0 %e jpt1 %e low %e up %e\n", ii, ibin, jpt0[ii], jpt1[ii], hxi->GetXaxis()->GetBinLowEdge(ibin), hxi->GetXaxis()->GetBinUpEdge(ibin)); exit(1);
    }
    htmp = hxi->ProjectionY(Form("hxi%d",ii), ibin, ibin); 
    //lout->Add(htmp); if(htmp->GetSumw2()==0) {htmp->Sumw2();} htmp->Scale(1./njet[ii]); 
    downscaleHist(htmp, njet[ii]);
    normHist(htmp);
    htmp->GetYaxis()->SetTitle("dN/d#xi/Njet");
    grtmp = XGLUtils::HistToGraph("grid3", htmp); lout->Add(grtmp);

    /*
    //should not be equal because hxi is filled per track, not per event!
    if( TMath::Abs(njet[ii] - hxi->Integral(ibin, ibin, 0, hxi->GetNbinsY()+1)) > EPSILON){
      printf("bad njet from hxi ii %d njet %e xi %e\n", ii, njet[ii], hxi->Integral(ibin, ibin, 0, hxi->GetNbinsY()+1)); exit(1);
    }
    */

    printf("ii %d ibin %d jpt0 %.1f jpt1 %.1f njet %e\n", ii, ibin, jpt0[ii], jpt1[ii], njet[ii]);

    TFile *fout=new TFile(Form("rawyield_TPC%d_JPT%.0f-%.0f.root", ktpc, jpt0[ii], jpt1[ii]),"recreate");
    lout->Write();
    fout->Save();
    fout->Close();
    delete fout;
  }
}

int main(int argc, char* argv[])
{
  for(Int_t ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }
  if(argc!=3){
    printf("argc!=3\n"); exit(1);
  }

  const Int_t tpccutmode = atoi(argv[2]);

  /* http://root.cern.ch/root/html/TH1.html
     When an histogram is created, a reference to it is automatically added to the list of in-memory objects for the current file or directory. This default behaviour can be changed by:
     h->SetDirectory(0);          for the current histogram h   
     TH1::AddDirectory(kFALSE);   sets a global switch disabling the reference  
     When the histogram is deleted, the reference to it is removed from the list of objects in memory. When a file is closed, all histograms in memory associated with this file are auto \
     matically deleted.    */
  //very important
  TH1::AddDirectory(kFALSE);

  TFile *_file0 = new TFile(argv[1]);
  printf("\nFile content:\n");
  _file0->ls();

  TString dn;
  TString ln;
  const Int_t tofid = 3;
  //mc
  //dn=Form("PWGJE_FragmentationFunction_clustersAOD_ANTIKT04_B0_Filter00272_Cut00150_Skip00_clustersAODMC2_ANTIKT04_B0_Filter00272_Cut00150_Skip00_AODMC2b_AODMCb_tpc%d_tof%d_cl0", tpccutmode, tofid);
  //ln=Form("fracfunc_clustersAOD_ANTIKT04_B0_Filter00272_Cut00150_Skip00_clustersAODMC2_ANTIKT04_B0_Filter00272_Cut00150_Skip00_AODMC2b_AODMCb_tpc%d_tof%d_cl0", tpccutmode, tofid);

  //GSI data
  //dn=Form("PWGJE_FragmentationFunction_jetsAOD_FASTJET04_B0_Filter00768_Cut00150_jetsAODMC2_FASTJET04_B0_Filter00768_Cut00150_AODMC2b_AODMCb_tpc%d_tof%d_cl0", tpccutmode, tofid);
  //ln=Form("fracfunc_jetsAOD_FASTJET04_B0_Filter00768_Cut00150_jetsAODMC2_FASTJET04_B0_Filter00768_Cut00150_AODMC2b_AODMCb_tpc%d_tof%d_cl0", tpccutmode, tofid);
  
  //LEGO
  dn=Form("PWGJE_FragmentationFunction_clustersAOD_ANTIKT04_B0_Filter00272_Cut00150_Skip00_noGenJets_trackTypeUndef_jetTypeUndef_tpc%d_tof%d_cl0", tpccutmode, tofid);
  ln=Form("fracfunc_clustersAOD_ANTIKT04_B0_Filter00272_Cut00150_Skip00_noGenJets_trackTypeUndef_jetTypeUndef_tpc%d_tof%d_cl0",tpccutmode, tofid);

  printf("please check\n%s\n%s\n", dn.Data(), ln.Data());

  TDirectory *fdir=_file0->GetDirectory(dn);
  if(!fdir){
    printf("no fdir!!\n");exit(1);
  }
  printf("\nDirectory content:\n");
  fdir->ls();

  TList * ll=(TList*)fdir->Get(ln);
  if(!ll){
    printf("no list!!\n");exit(1);
  }
  printf("\nList content:\n");
  ll->ls();

  getRawYield(ll, tpccutmode);

  printf("\ngetRawYield done\n");
  return 0;
}
